#include "config.h"

#include "lstatfit.h"

#define GAP	0.00001

ecm_double
exp_goodness(msl_fit_t *fit)
{
	ecm_uint32	i;
	ecm_double	*yvalues = fit->yvalues;
	ecm_double	sq_diffs = 0;
	ecm_double	lambda;

	lambda = gsl_vector_get(fit->s->x, 0);

	for (i = 0; i < fit->n; i++) {
		ecm_double	prob;

		prob = 1 - exp(-lambda * i);

		sq_diffs += pow(yvalues[i] - prob, 2);
	}

	return sq_diffs;
}

void
exp_dist_report(FILE *fp, msl_fit_t *fit, ecm_uint32 xdiff)
{
	ecm_uint32	i;
	ecm_double	*xvalues = fit->xvalues;
	ecm_double	lambda;

	lambda = gsl_vector_get(fit->s->x, 0);

	ecm_fprintf(fp, "# lambda:%lf\n", lambda);
	
	for (i = 0; i < fit->n; i++) {
		double	t = xvalues[i];
		ecm_double	prob;

		prob = 1 - exp(-lambda * t);
		ecm_fprintf(fp, "%8u %lf\n", t + xdiff, prob);
	}
}

int
exp_dist_f(const gsl_vector *x, void *data, gsl_vector *f)
{
	msl_fit_t	*fit = (msl_fit_t *)data;
	ecm_uint32	n = fit->n;
	ecm_double	*xvalues = fit->xvalues;
	ecm_double	*yvalues = fit->yvalues;
     
	double	lambda = gsl_vector_get(x, 0);

	ecm_uint32	i;
     
	for (i = 0; i < n; i++)	{
		/* Model Yi = lambda * exp(-lambda * i) */
		/* Model Yi = 1 - exp(-lambda * i) */
		double	t = xvalues[i];
		double	Yi = 1 - exp(-lambda * t);

		gsl_vector_set(f, i, Yi - yvalues[i]);
	}
     
	return GSL_SUCCESS;
}

int
exp_dist_df(const gsl_vector *x, void *data, gsl_matrix *J)
{
	msl_fit_t	*fit = (msl_fit_t *)data;
	ecm_uint32	n = fit->n;
	ecm_double	*xvalues = fit->xvalues;

	double	lambda = gsl_vector_get(x, 0);

	size_t	i;
     
	for (i = 0; i < n; i++)	{
		double	t = xvalues[i];
		double	e = exp(-lambda * t);
		double	dlambda;

		dlambda = t * e;
		
		gsl_matrix_set(J, i, 0, dlambda);
	}
	return GSL_SUCCESS;
}

